Appendix 2: Landscapes Landuse and forest cover at local and landscape scale
load(here("data","datasets.Rdata")) # datasets
load("models_outputs/models_local_cover.Rdata") # models outputs
load("models_outputs/models_landscape_cover.Rdata") # models outputs
source("scripts/r2_auxiliary_function.R")# variance inflation factor function
vif.mer <- function (fit) {
## adapted from rms::vif
v <- vcov(fit)
nam <- names(fixef(fit))
## exclude intercepts
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
v
}This appendix is a description of the landuse and matrix of the landscapes, togethe with a more specific description of forest cover variables at both local and landscape scales. We present the baseline models for the selection of the best scale for the local forest cover variable for each dataset. For the local scale, we measured the percentage of forest cover within buffers of 400, 600 and 800 m around each sampling site. For the landscape forest cover, we used the 2 km buffer around the landscape centroid.
1. Landuse and matrix description in landscapes
lu <- landuse %>%
mutate(landscape = fct_reorder(landscape,forest,max)) %>%
pivot_longer(3:10, names_to = "landuse", values_to = "value") %>%
mutate(landuse = fct_relevel(landuse,
"forest", "pasture", "coffee", "eucalyptus",
"sugarcane", "urban", "water","other")) Percentage of each landuse type per landscape
# label colors for each landuse
cores <- c("#549E3E", "#FDF5B5", "#D97853", # for pas coffee
"#BCDE93", "#FFBBFF", "#BDB39A", # euc sugar water
"#7AC5CD", "#8B7D6B") # water other
lu %>% mutate(landuse = fct_rev(landuse)) %>%
ggplot(aes(y=landscape, x=value, fill=landuse)) +
geom_col() +
facet_grid(matrix~., scales="free_y",space='free_y', switch="y",
labeller = as_labeller(c("high-quality" = "high-quality landscapes",
"low-quality" = "Low-quality landscapes"))) +
scale_fill_manual(values=rev(cores)) +
theme(legend.position = "bottom",
strip.placement= "outside")+
xlab("Landuse (%)") + ylab("") +
ggtitle("Landuse composition")Proportions of landuse types among high- and low-quality landscapes
ggplot(lu, aes(x=matrix, y=value,fill=landuse, col=landuse)) +
facet_wrap(~landuse, ncol=4) +
geom_violin(fill="white")+
geom_dotplot(stackdir = "center", dotsize=1.5,
binaxis = "y") +
stat_summary(fun=median, size=0.1, fatten=5, col="black",
geom="crossbar", width=0.5)+
scale_fill_manual(values=cores )+
scale_color_manual(values=cores ) +
theme(legend.position = "none") +
xlab("Matrix quality") + ylab("Landuse (%)")Summary table
lu %>% group_by(landuse, matrix) %>% summarise(min = min(value),
mean=mean(value),
sd = sd(value), median=median(value),
max = max(value)) %>%
kable(digits = 1)| landuse | matrix | min | mean | sd | median | max |
|---|---|---|---|---|---|---|
| forest | high-quality | 12.0 | 30.8 | 12.4 | 29.0 | 55.5 |
| forest | low-quality | 6.9 | 25.3 | 11.9 | 26.8 | 46.0 |
| pasture | high-quality | 20.4 | 35.5 | 10.2 | 34.7 | 49.5 |
| pasture | low-quality | 31.7 | 46.9 | 11.1 | 44.8 | 69.1 |
| coffee | high-quality | 12.4 | 21.7 | 6.8 | 19.5 | 32.9 |
| coffee | low-quality | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| eucalyptus | high-quality | 0.0 | 1.6 | 3.2 | 0.4 | 10.1 |
| eucalyptus | low-quality | 4.3 | 18.1 | 8.7 | 16.9 | 37.7 |
| sugarcane | high-quality | 0.0 | 4.2 | 7.9 | 0.0 | 20.6 |
| sugarcane | low-quality | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| urban | high-quality | 0.8 | 2.1 | 1.7 | 1.4 | 5.5 |
| urban | low-quality | 0.7 | 7.4 | 7.5 | 3.8 | 26.6 |
| water | high-quality | 0.0 | 2.5 | 4.6 | 0.3 | 12.7 |
| water | low-quality | 0.0 | 0.9 | 1.9 | 0.1 | 6.2 |
| other | high-quality | 0.0 | 1.6 | 1.1 | 1.5 | 3.4 |
| other | low-quality | 0.0 | 1.4 | 2.5 | 0.0 | 8.0 |
Principal Coordinate Analysis
To compare landuse matrices in both high- and low-quality landscapes, we performed a Principal Coordiante analysis with the proportions of each landuse type. - PCoA is more advisable given the nature of the data (sum up 100%, compositional data), we used gower distance tranformation to create a distance matrix among landscapes
cor <- as.numeric(as.factor(landuse$matrix)) # landuse colors
pcoa <- capscale (landuse[,3:10]~1, dist= "gower", scaling = 1)
plot(pcoa, main = 'PCoA (MDS)', display="si", scaling=2, type="points",
ylim=c(-1.1,1),
col=cor, pch=16)
text (pcoa, display = 'sp', col="blue", scaling=3)
points (pcoa, display = 'si', col=cor, pch=16,scaling=2, cex=1.2)Landuse heterogeneity - diversity index
Comparing matrix
landuse$simp.div <- diversity(landuse[,3:10], index='simpson')
landuse$shan.div <- diversity(landuse[,3:10], index='shannon')
ggplot(landuse, aes(x=matrix, y=simp.div)) + geom_boxplot()+
geom_dotplot(stackdir = "center", dotsize=1,binaxis = "y") +
ggplot(landuse, aes(x=matrix, y=shan.div)) + geom_boxplot() +
geom_dotplot(stackdir = "center", dotsize=1, binaxis = "y")Testing differences in diversity indexes between matrix quality
summary(lm(simp.div ~matrix, data=landuse))##
## Call:
## lm(formula = simp.div ~ matrix, data = landuse)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.13997 -0.04294 0.01388 0.04218 0.08816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.69288 0.01761 39.353 <2e-16 ***
## matrixlow-quality -0.05281 0.02342 -2.255 0.0349 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05568 on 21 degrees of freedom
## Multiple R-squared: 0.1949, Adjusted R-squared: 0.1566
## F-statistic: 5.085 on 1 and 21 DF, p-value: 0.03494
summary(lm(shan.div ~matrix, data=landuse))##
## Call:
## lm(formula = shan.div ~ matrix, data = landuse)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19609 -0.10582 -0.01792 0.09598 0.34078
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.34359 0.04329 31.036 <2e-16 ***
## matrixlow-quality -0.15295 0.05758 -2.656 0.0148 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1369 on 21 degrees of freedom
## Multiple R-squared: 0.2515, Adjusted R-squared: 0.2158
## F-statistic: 7.055 on 1 and 21 DF, p-value: 0.01478
Correlation among proportions of landuse types
ggpairs(landuse, columns = 3:7, ggplot2::aes(colour=matrix),
lower = list(continuous = "smooth"))2. Relationships among forest cover variables
We calculated Pearson correlation coefficients for forest cover variables in each matrix quality region (Figure S2.). Also, we ploted the range of local forest cover (400 m) within the landscapes to see how local forest cover varies among landscapes in both regions (Figure S2.).
high.cor <- env[which(env$matrix == "high_quality"), c(10,11,12,9)]
colnames(high.cor) <- c("Local 400m", "Local 600m",
"Local 800m","Landscape 2km")
low.cor <- env[which(env$matrix == "low_quality"), c(10,11,12,9)]
colnames(low.cor) <- c("Local 400m", "Local 600m",
"Local 800m","Landscape 2km")
par(mfrow=c(1,2))
corrplot(cor(high.cor),method="number",type="lower", diag=F,
order="original", cl.pos = "n",
mar = c(0,0,0,0))
mtext("High-quality matrix", cex=1.2)
corrplot(cor(low.cor),method="number",type="lower", diag=F,
order="original", cl.pos = "n")
mtext("Low-quality matrix", cex=1.2)Correlations among forest cover variables in the high-quality (left) and low-quality matrix (right) landscapes.
ggplot(env, aes(x=forest_site400, y=forest_land)) +
facet_grid(~matrix, labeller = as_labeller(c(high_quality = "High-quality matrix", low_quality = "Low-quality matrix"))) +
geom_point() +
geom_line(aes(shape=landscape)) +
xlab("% local forest cover (400 m)") +
ylab("% landscape forest cover (2 km)") +
ylim(5,60)Range of local forest cover (buffer 400 m) within each landscape (buffer 2 km) for both landscapes with different matrix quality. Each line represent a landscape and the dots area the local forest cover for each sampling site.
3. Scale of effects for local forest cover
We ran different models with each local forest cover variable and selected the scale of effect using AIC model selection and the R2 of the models. The models follow the specification presented in the paper (Modeling section), except that here we did not include trait variables, i.e. we only modeled the occurrence of species according to forest cover.
We used lme4 package (Bates et al. (2015)) to perform a GLMM with binomial (proportion) distribution. An example of the code for each assemblage is as follows:
model <- glmer(cbind(occor, n.visit-occor) ~
local.cover + (local.cover|sp) +
(1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.spe)
mhigh.spe4a <- glmer(cbind(occor, n.visit-occor) ~
forest_site400 +
(forest_site400 |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mhigh.spe6a <- glmer(cbind(occor, n.visit-occor) ~
forest_site600 +
(forest_site600 |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mhigh.spe8a <- glmer(cbind(occor, n.visit-occor) ~
forest_site800 +
(forest_site800|sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mlow.spe4a <- glmer(cbind(occor, n.visit-occor) ~
forest_site400 +
(forest_site400 |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=low.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mlow.spe6a <- glmer(cbind(occor, n.visit-occor) ~
forest_site600 +
(forest_site600 |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=low.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mlow.spe8a <- glmer(cbind(occor, n.visit-occor) ~
forest_site800 +
(forest_site800|sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=low.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mhigh.gen4a <- glmer(cbind(occor, n.visit-occor) ~
forest_site400 +
(forest_site400 |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.gen,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mhigh.gen6a <- glmer(cbind(occor, n.visit-occor) ~
forest_site600 +
(forest_site600 |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.gen,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mhigh.gen8a <- glmer(cbind(occor, n.visit-occor) ~
forest_site800 +
(forest_site800|sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.gen,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mlow.gen4a <- glmer(cbind(occor, n.visit-occor) ~
forest_site400 +
(forest_site400 |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=low.gen,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mlow.gen6a <- glmer(cbind(occor, n.visit-occor) ~
forest_site600 +
(forest_site600 |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=low.gen,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mlow.gen8a <- glmer(cbind(occor, n.visit-occor) ~
forest_site800 +
(forest_site800|sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=low.gen,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
save(mhigh.spe4a,mhigh.spe6a,mhigh.spe8a,
mhigh.gen4a, mhigh.gen6a, mhigh.gen8a,
mlow.spe4a, mlow.spe6a, mlow.spe8a,
mlow.gen4a,mlow.gen6a,mlow.gen8a,
file="models_outputs/models_local_cover.Rdata")taic.ce <- AICtab(mhigh.spe4a, mhigh.spe6a, mhigh.spe8a, sort=F)
class(taic.ce) <- "data.frame"
taic.pe <- AICtab(mlow.spe4a, mlow.spe6a, mlow.spe8a, sort=F)
class(taic.pe) <- "data.frame"
taic.cg <- AICtab(mhigh.gen4a, mhigh.gen6a, mhigh.gen8a, sort=F)
class(taic.cg) <- "data.frame"
taic.pg <- AICtab(mlow.gen4a, mlow.gen6a, mlow.gen8a, sort=F)
class(taic.pg) <- "data.frame"
taic <- rbind(taic.pe, taic.ce, taic.pg, taic.cg)taic$modelo <- rep(c("400m", "600m", "800m"), 4)
taic$dataset <- rep(c("low.spe", "high.spe", "low.gen", "high.gen"), each=3)
ttaic <- taic %>% dplyr::select(dataset,modelo,dAIC) %>% spread(dataset, dAIC) %>% select(modelo, high.spe, low.spe, high.gen, low.gen)
colnames(ttaic) <- c("model", "Low-quality", "High-quality", "Low-quality", "High-quality")r2high.spe <- data.frame(r.quad(mhigh.spe4a, "forest_site400"),
m6a = r.quad(mhigh.spe6a, c("forest_site600"))[,2],
m8a = r.quad(mhigh.spe8a, "forest_site800")[,2]) %>%
mutate_if(is.numeric, round, digits=3)
colnames(r2high.spe)[2] <- "m4a"
rownames(r2high.spe) <- r2high.spe$componente
r2low.spe <- data.frame(r.quad(mlow.spe4a, "forest_site400"),
m6a = r.quad(mlow.spe6a, c("forest_site600"))[,2],
m8a = r.quad(mlow.spe8a, "forest_site800")[,2]) %>%
mutate_if(is.numeric, round, digits=3)
colnames(r2low.spe)[2] <- "m4a"
rownames(r2low.spe) <- r2low.spe$componente
r2high.gen <- data.frame(r.quad(mhigh.gen4a, "forest_site400"),
m6a = r.quad(mhigh.gen6a, c("forest_site600"))[,2],
m8a = r.quad(mhigh.gen8a, "forest_site800")[,2]) %>%
mutate_if(is.numeric, round, digits=3)
colnames(r2high.gen)[2] <- "m4a"
rownames(r2high.gen) <- r2high.gen$componente
r2low.gen <- data.frame(r.quad(mlow.gen4a, "forest_site400"),
m6a = r.quad(mlow.gen6a, c("forest_site600"))[,2],
m8a = r.quad(mlow.gen8a, "forest_site800")[,2]) %>%
mutate_if(is.numeric, round, digits=3)
colnames(r2low.gen)[2] <- "m4a"
rownames(r2low.gen) <- r2low.gen$componente
tudo <- rbind( t(r2low.spe[,-1]),t(r2high.spe[,-1]), t(r2low.gen[,-1]), t(r2high.gen[,-1]))
tudo <- tudo*100tabis <- cbind(tudo,taic) %>% mutate(dAIC = round(dAIC,2))
colnames(tabis) <- c("Total", "fixed","env.sp", "lands.sp", "site.sp", "lands", "site", "dAIC", "df", "Model","dataset")
kable(tabis[,c(10,1:9)], "latex", booktabs=T, row.names = F,
caption= "Overal and marginal r-squares and model comparisons with Akaike Information Criterion (AIC) for models with different local forest cover scales as predictor for the specialist and generalists species in both regions with different matrix qualities. For the terms see Table 1 (main text). dAIC is the difference in Akaike Information Criterion to the best model; df are the degrees of freedom.") %>%
kable_styling() %>%
group_rows("Forest specialist species", 1, 6) %>%
group_rows(" Low-quality matrix", 1,3) %>%
group_rows(" High-quality matrix", 4,6) %>%
group_rows("Forest generalist species", 7, 12) %>%
group_rows(" Low-quality matrix", 7,9) %>%
group_rows(" High-quality matrix", 10,12) %>%
add_header_above(c(" " = 1, " " = 7, "AIC"= 2))pforest_site400 <- data.frame(forest_site400 = seq(min(high.spe$forest_site400), max(high.spe$forest_site400), length.out=20))
pforest_site400$pred <- predict(mhigh.spe4a, newdata=pforest_site400, re.form=NA, type="response")
pforest_site400$cob <- pforest_site400$forest_site400*sd(high.spe$forest_site400Orig) + mean(high.spe$forest_site400Orig)
pforest_site600 <- data.frame(forest_site600 = seq(min(high.spe$forest_site600), max(high.spe$forest_site600), length.out=20))
pforest_site600$pred <- predict(mhigh.spe6a, newdata=pforest_site600, re.form=NA, type="response")
pforest_site600$cob <- pforest_site600$forest_site600*sd(high.spe$forest_site600Orig) + mean(high.spe$forest_site600Orig)
pforest_site800 <- data.frame(forest_site800 = seq(min(high.spe$forest_site800), max(high.spe$forest_site800), length.out=20))
pforest_site800$pred <- predict(mhigh.spe8a, newdata=pforest_site800, re.form=NA, type="response")
pforest_site800$cob<- pforest_site800$forest_site800*sd(high.spe$forest_site800Orig) + mean(high.spe$forest_site800Orig)
fhigh.spe <- rbind(pforest_site400[,2:3], pforest_site600[,2:3], pforest_site800[,2:3]) %>%
mutate(modelo = rep(c("forest_site400","forest_site600","forest_site800"), each=20),
dataset = "high.spe")
#ggplot(fig, aes(x=cob, y=pred, col=modelo)) + geom_line()pforest_site400 <- data.frame(forest_site400 = seq(min(low.spe$forest_site400), max(low.spe$forest_site400), length.out=20))
pforest_site400$pred <- predict(mlow.spe4a, newdata=pforest_site400, re.form=NA, type="response")
pforest_site400$cob <- pforest_site400$forest_site400*sd(low.spe$forest_site400Orig) + mean(low.spe$forest_site400Orig)
pforest_site600 <- data.frame(forest_site600 = seq(min(low.spe$forest_site600), max(low.spe$forest_site600), length.out=20))
pforest_site600$pred <- predict(mlow.spe6a, newdata=pforest_site600, re.form=NA, type="response")
pforest_site600$cob <- pforest_site600$forest_site600*sd(low.spe$forest_site600Orig) + mean(low.spe$forest_site600Orig)
pforest_site800 <- data.frame(forest_site800 = seq(min(low.spe$forest_site800), max(low.spe$forest_site800), length.out=20))
pforest_site800$pred <- predict(mlow.spe8a, newdata=pforest_site800, re.form=NA, type="response")
pforest_site800$cob<- pforest_site800$forest_site800*sd(low.spe$forest_site800Orig) + mean(low.spe$forest_site800Orig)
flow.spe <- rbind(pforest_site400[,2:3], pforest_site600[,2:3], pforest_site800[,2:3]) %>%
mutate(modelo = rep(c("forest_site400","forest_site600", "forest_site800"), each=20),
dataset = "low.spe")pforest_site400 <- data.frame(forest_site400 = seq(min(high.gen$forest_site400), max(high.gen$forest_site400), length.out=20))
pforest_site400$pred <- predict(mhigh.gen4a, newdata=pforest_site400, re.form=NA, type="response")
pforest_site400$cob <- pforest_site400$forest_site400*sd(high.gen$forest_site400Orig) + mean(high.gen$forest_site400Orig)
pforest_site600 <- data.frame(forest_site600 = seq(min(high.gen$forest_site600), max(high.gen$forest_site600), length.out=20))
pforest_site600$pred <- predict(mhigh.gen6a, newdata=pforest_site600, re.form=NA, type="response")
pforest_site600$cob <- pforest_site600$forest_site600*sd(high.gen$forest_site600Orig) + mean(high.gen$forest_site600Orig)
pforest_site800 <- data.frame(forest_site800 = seq(min(high.gen$forest_site800), max(high.gen$forest_site800), length.out=20))
pforest_site800$pred <- predict(mhigh.gen8a, newdata=pforest_site800, re.form=NA, type="response")
pforest_site800$cob<- pforest_site800$forest_site800*sd(high.gen$forest_site800Orig) + mean(high.gen$forest_site800Orig)
fhigh.gen <- rbind(pforest_site400[,2:3], pforest_site600[,2:3], pforest_site800[,2:3]) %>%
mutate(modelo = rep(c("forest_site400","forest_site600","forest_site800"), each=20),
dataset = "high.gen")pforest_site400 <- data.frame(forest_site400 = seq(min(low.gen$forest_site400), max(low.gen$forest_site400), length.out=20))
pforest_site400$pred <- predict(mlow.gen4a, newdata=pforest_site400, re.form=NA, type="response")
pforest_site400$cob <- pforest_site400$forest_site400*sd(low.gen$forest_site400Orig) + mean(low.gen$forest_site400Orig)
pforest_site600 <- data.frame(forest_site600 = seq(min(low.gen$forest_site600), max(low.gen$forest_site600), length.out=20))
pforest_site600$pred <- predict(mlow.gen6a, newdata=pforest_site600, re.form=NA, type="response")
pforest_site600$cob <- pforest_site600$forest_site600*sd(low.gen$forest_site600Orig) + mean(low.gen$forest_site600Orig)
pforest_site800 <- data.frame(forest_site800 = seq(min(low.gen$forest_site800), max(low.gen$forest_site800), length.out=20))
pforest_site800$pred <- predict(mlow.gen8a, newdata=pforest_site800, re.form=NA, type="response")
pforest_site800$cob<- pforest_site800$forest_site800*sd(low.gen$forest_site800Orig) + mean(low.gen$forest_site800Orig)
flow.gen <- rbind(pforest_site400[,2:3], pforest_site600[,2:3], pforest_site800[,2:3]) %>%
mutate(modelo = rep(c("forest_site400", "forest_site600", "forest_site800"),each=20),
dataset="low.gen")In Figure S2., we present the occurrence probability predicted for the models with different local forest cover scales for all the assemblages. Predictions were quite similar and decreased with forest cover for the specialists, especially in the low-quality matrix region, and increased or remained flat for the generalists.
labs <- c(high = "High-quality matrix", low="Low-quality matrix",
spe = "Specialists", gen = "Generalists")
rbind(fhigh.spe, flow.spe,fhigh.gen, flow.gen) %>%
separate(dataset, c("matrix", "habitat")) %>%
ggplot(aes(x=cob, y=pred, col=modelo)) + geom_line()+
facet_grid(habitat~matrix, labeller=as_labeller(labs), scales = "free") +
scale_color_discrete("Scale", labels=c("400 m", "600 m", "800 m")) +
scale_x_continuous(limits = c(0,80), name="Local forest cover (%)")+
scale_y_continuous("Occurrence probability",
breaks=c(0,0.02,0.04,0.06,0.08,0.1),
expand=expand_scale(add=0), limits=c(0,0.1)) +
theme(panel.spacing = unit(5, "mm")) +
geom_hline(yintercept = 0, size=0.8)Predictions of the models with different local forest cover scales (lines) for specialists and generalists in both regions.
Residual correlations among species
We evaluated the residuals by Kendall correlations among species and among sites for the 400 m models using the predictions for site:sp random effects (Observation Level Random Effect), following the code provided by Miller, Damschen & Ives (2018). Codes for the species names are presented in the dataset available.
Below we show models residual correlation plots for the specialists in the high-quality matrix landscape. All the other assemblages presented similar results.
nsites = 40
nspp = length(unique(high.spe$sp))
dat <- high.spe
dat$resid <- as.matrix(ranef(mhigh.spe4a)$`site:sp`)
X <- matrix(dat$resid, nrow=nsites, ncol=nspp, byrow=F)
rownames(X) <- unique(dat$site)
colnames(X) <- unique(dat$sp)
# species correlations
corrs.sp <- cor(X, method="kendall")
for(i in 1:nspp) corrs.sp[i,i] <- NA
# site correlations
corrs.site <- cor(t(X), method="kendall")
for(i in 1:nsites) corrs.site[i,i] <- NARange of species correlations: -0.41, 0.46.
Range of sites correlations: -0.25, 0.31.
corrplot(corrs.sp, method="square", type="upper",diag=F,tl.cex=0.5,cl.cex=0.5)Species residual Kendall correlations for the specialist species in the coffee matrix region.
corrplot(corrs.site, method="square", type="upper",diag=F,tl.cex=0.5,cl.cex=0.5)Sites residual Kendall correlations for the specialist species in the coffee matrix region.
x.sp <- matrix(corrs.sp, ncol=1)
x.site <- matrix(corrs.site, ncol=1)
x.sp <- x.sp[!is.na(x.sp)]
x.site <- x.site[!is.na(x.site)]
# Figure histogram
par(mfrow=c(1,2), mai=c(1,1,.5,.2))
hist(corrs.sp, main="", probability=T, breaks=.1*(-10:10), xlab="Correlation coefficients", ylab="Probability", ylim=c(0,4), cex.lab=1.2)
mtext("Species", side=3,cex=1.2)
hist(corrs.site, main="", probability=T, breaks=.1*(-10:10), xlab="Correlation coefficients", ylab="Probability", ylim=c(0,4), cex.lab=1.2)
mtext("Sites", side=3, cex=1.2)Histograms of the residual Kendall correlations for the specialists species in the coffee matrix region.
4. Including landscape forest cover
After selecting the local forest cover of 400 m radius buffer around each site for all datasets, we included the landscape forest cover (2 km radius buffer around the centroid of the landscape) in the model.
The R syntax example of this model area as follows:
model <- glmer(cbind(occor, n.visit-occor) ~
local.400 + landscape.2k +
(local.400 + landscape.2k |sp) +
(1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.spe)
mhigh.spe42 <- glmer(cbind(occor, n.visit-occor) ~
forest_site400 + forest_land +
(forest_site400 + forest_land |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mlow.spe42 <- glmer(cbind(occor, n.visit-occor) ~
forest_site400 + forest_land +
(forest_site400 + forest_land |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=low.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mhigh.gen42 <- glmer(cbind(occor, n.visit-occor) ~
forest_site400 + forest_land +
(forest_site400 + forest_land |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=high.gen,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
mlow.gen42 <- glmer(cbind(occor, n.visit-occor) ~
forest_site400 + forest_land +
(forest_site400 + forest_land |sp) + (1|landscape:sp) + (1|site:sp) +
(1|landscape) + (1|site),
family=binomial, data=low.gen,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
save(mhigh.spe42, mlow.spe42, mhigh.gen42, mlow.gen42,
file="models_outputs/models_landscape_cover.Rdata")Before analysing results, we evaluated possible colinearity between local and landscape forest cover using the Variance Inflation Factor with the code provided by John Lefcheck (https://jonlefcheck.net/2012/12/28/dealing-with-multicollinearity-using-variance-inflation-factors/). With VIF we found no evidence of collinearity between the forest cover scales (Table S2.).
varif <- rbind(vif.mer(mhigh.spe42),vif.mer(mlow.spe42),vif.mer(mhigh.gen42),
vif.mer(mlow.gen42))
rownames(varif) <- c("Coffee", "Pasture", "Coffe", "Pasture")
colnames(varif) <- c("Local", "Landscape")
kable(varif, "latex", booktabs=T, caption= "Variance Inflation Factor index for the variables of local forest cover and landscape forest cover.", digits=2) %>%
group_rows("Specialists", 1, 2) %>%
group_rows("Generalists", 3,4) Predictions of the models are present in Figure S2.. It is important to notice the differences in 20 and 40% landscape forest cover predictions for the specialists in the low-quality matrix.
pforest_site400 <- data.frame(forest_site400 = seq(min(high.spe$forest_site400),
max(high.spe$forest_site400), length.out=20))
pforest_site400$pred <- predict(mhigh.spe4a, newdata=pforest_site400, re.form=NA, type="response")
pforest_site400$cob <- pforest_site400$forest_site400*sd(high.spe$forest_site400Orig) + mean(high.spe$forest_site400Orig)
porc <- c(20, 40)
p42 <- expand.grid(forest_site400 = seq(min(high.spe$forest_site400), max(high.spe$forest_site400), length.out=10),
forest_land = (porc - mean(high.spe$forest_landOrig))/sd(high.spe$forest_landOrig))
p42$pred <- predict(mhigh.spe42, newdata=p42, re.form=NA, type="response")
p42$cob <- p42$forest_site400*sd(high.spe$forest_site400Orig) + mean(high.spe$forest_site400Orig)
p42$cob.paisa <- p42$forest_land*sd(high.spe$forest_landOrig) + mean(high.spe$forest_landOrig)
loc <- pforest_site400[,2:3] %>%
mutate(cob.paisa=NA, modelo = "4a")
pais <- p42[,3:5] %>%
mutate(modelo ="42")
fhigh.spe <- rbind(loc, pais) %>% mutate(dataset="high.spe")pforest_site400 <- data.frame(forest_site400 = seq(min(low.spe$forest_site400), max(low.spe$forest_site400), length.out=20))
pforest_site400$pred <- predict(mlow.spe4a, newdata=pforest_site400, re.form=NA, type="response")
pforest_site400$cob <- pforest_site400$forest_site400*sd(low.spe$forest_site400Orig) + mean(low.spe$forest_site400Orig)
porc <- c(20, 40)
p42 <- expand.grid(forest_site400 = seq(min(low.spe$forest_site400), max(low.spe$forest_site400), length.out=10),
forest_land = (porc - mean(low.spe$forest_landOrig))/sd(low.spe$forest_landOrig))
p42$pred <- predict(mlow.spe42, newdata=p42, re.form=NA, type="response")
p42$cob <- p42$forest_site400*sd(low.spe$forest_site400Orig) + mean(low.spe$forest_site400Orig)
p42$cob.paisa <- p42$forest_land*sd(low.spe$forest_landOrig) + mean(low.spe$forest_landOrig)
loc <- pforest_site400[,2:3] %>%
mutate(cob.paisa=NA, modelo = "4a")
pais <- p42[,3:5] %>%
mutate(modelo ="42")
flow.spe <- rbind(loc, pais) %>% mutate(dataset="low.spe")pforest_site400 <- data.frame(forest_site400 = seq(min(high.gen$forest_site400), max(high.gen$forest_site400), length.out=20))
pforest_site400$pred <- predict(mhigh.gen4a, newdata=pforest_site400, re.form=NA, type="response")
pforest_site400$cob <- pforest_site400$forest_site400*sd(high.gen$forest_site400Orig) + mean(high.gen$forest_site400Orig)
porc <- c(20, 40)
porc <- c(20, 40)
p42 <- expand.grid(forest_site400 = seq(min(high.gen$forest_site400), max(high.gen$forest_site400), length.out=10),
forest_land = (porc - mean(high.gen$forest_landOrig))/sd(high.gen$forest_landOrig))
p42$pred <- predict(mhigh.gen42, newdata=p42, re.form=NA, type="response")
p42$cob <- p42$forest_site400*sd(high.gen$forest_site400Orig) + mean(high.gen$forest_site400Orig)
p42$cob.paisa <- p42$forest_land*sd(high.gen$forest_landOrig) + mean(high.gen$forest_landOrig)
loc <- pforest_site400[,2:3] %>%
mutate(cob.paisa=NA, modelo = "4a")
pais <- p42[,3:5] %>%
mutate(modelo ="42")
fhigh.gen <- rbind(loc, pais) %>% mutate(dataset="high.gen")pforest_site400 <- data.frame(forest_site400 = seq(min(low.gen$forest_site400), max(low.gen$forest_site400), length.out=20))
pforest_site400$pred <- predict(mlow.gen4a, newdata=pforest_site400, re.form=NA, type="response")
pforest_site400$cob <- pforest_site400$forest_site400*sd(low.gen$forest_site400Orig) + mean(low.gen$forest_site400Orig)
porc <- c(20, 40)
p42 <- expand.grid(forest_site400 = seq(min(low.gen$forest_site400), max(low.gen$forest_site400), length.out=10),
forest_land = (porc - mean(low.gen$forest_landOrig))/sd(low.gen$forest_landOrig))
p42$pred <- predict(mlow.gen42, newdata=p42, re.form=NA, type="response")
p42$cob <- p42$forest_site400*sd(low.gen$forest_site400Orig) + mean(low.gen$forest_site400Orig)
p42$cob.paisa <- p42$forest_land*sd(low.gen$forest_landOrig) + mean(low.gen$forest_landOrig)
loc <- pforest_site400[,2:3] %>%
mutate(cob.paisa=NA, modelo = "4a")
pais <- p42[,3:5] %>%
mutate(modelo ="42")
flow.gen <- rbind(loc, pais) %>% mutate(dataset="low.gen")labs <- c(high = "High-quality matrix", low = "Low-quality matrix",
spe = "Specialists", gen = "Generalists")
rbind(fhigh.spe, flow.spe, fhigh.gen, flow.gen) %>%
separate(dataset, c("matrix", "habitat")) %>%
ggplot(aes(x=cob, y=pred, col=as.factor(cob.paisa))) +
geom_line() +
facet_grid(habitat~matrix, labeller=as_labeller(labs)) +
scale_color_discrete("Landscape cover", labels=c("20%", "40%", "only local model")) +
scale_x_continuous(limits = c(0,80), name="Local forest cover 400m (%)") +
scale_y_continuous("Occurrence probability",
expand=expand_scale(add=0), limits=c(0,0.16)) +
theme(panel.spacing = unit(5, "mm")) +
geom_hline(yintercept = 0, size=0.8)Predictions of the models without (gray lines) and with landscape forest cover scales (20 percent cover in red and 40 percent cover in blue lines) for specialists and generalists in both regions.